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A strongly coupled quark-gluon plasma (QGP) of heavy constituent quasi-particles is studied 
by a path-integral Monte-Carlo method. This approach is a quantum generalization of the model 
developed by Gelman, Shuryak and Zahed. It is shown that this method is able to reproduce the 
QCD lattice equation of state and also yields valuable insight into the internal structure of the QGP. 
The results indicate that the QGP reveals liquid-like rather than gas-like properties. At temperatures 
just above the critical one it was found that bound quark-antiquark states still survive. These states 
are bound by effective string-like forces and turns out to be colorless. At the temperature as large 
as twice the critical one no bound states are observed. Quantum effects turned out to be of prime 
importance in these simulations. 



I. INTRODUCTION 



Contributions of A.B. Migdal to physics of the last century are difRcult to overestimate. These include the Landau- 
Migdal-Pomeranchuk effect, theory of finite Fermi systems, pionic degrees of freedom in nuclei, physics of strong 
fields, etc. This results entered many textbooks and found numerous applications. A.B. Migdal was one of the first 
who started to think about possible existence of an anomalous nuclear matter Nowadays this concept is, in 
particular, associated with the phase transition into the QGP. While the possibility of pion condensation is still 
vividly discussed in astrophysical applications, see, e.g., 0. 

Investigation of properties of the QGP is one of the main challenges of strong-interaction physics, both theoretically 
and experimentally. Many features of this matter were experimentally discovered at the Relativistic Heavy Ion Collider 
(RHIC) at Brookhaven. The most striking result, obtained from analysis of these experimental data is that the 
deconfined quark-gluon matter behaves as an almost perfect fluid rather than a perfect gas, as it could be expected 
from the asymptotic freedom. 

There are various approaches to studying QGP. Each approach has its advantages and disadvantages. The most 
fundamental way to compute properties of the strongly interacting matter is provided by the lattice QCD 0- 
0. Interpretation of these very complicated computations requires application of various QCD motivated, albeit 
schematic, models simulating various aspects of the full theory. Moreover, suclr models are needed in cases when the 
lattice QCD fails, e.g. at large baryon chemical potentials and out of equilibrium. While some progress has been 
achieved in the recent years, we are still far away from having a satisfactory understanding of the QlGP dynamics. 

A semi-classical approximation, based on a point like quasi-particle picture has been introduced in [7|. It is expected 
that the main features of non-Abelian plasmas can be understood in simple semi-classical terms without the difficulties 
inherent to a full quantum field theoretical analysis. Independently the same ideas were implemented in terms of 
molecular dynamics (MD) Q. Recently this MD approach was further developed in a series of works The MD 

allowed one to treat soft processes in the QGP which are not accessible by perturbative means. 

A strongly correlated behavior of the QGP is expected to show up in long-ranged spatial correlations of quarks 
and gluons which, in fact, may give rise to liquid-like and, possibly, solid-like structures. This expectation is based on 
a very similar behavior observed in electrodynamic plasmas [ill . [T3 . This similarity was exploited to formulate a 
classical non-relativistic model of a color Coulomb interacting QGP [0| which was numerically analyzed by classical 
MD simulations. Quantum effects were either neglected or included phenomenologically via a short-range repulsive 
correction to the pair potential. Such a rough model may become a critical issue at high densities, where quantum 
effects strongly affects properties of the QGP. Similar models have been used in electrodynamic plasmas and showed 
poor behavior in the region of strong wave function overlap, in particular at the Mott density. For temperatures and 
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densities of the QGP considered in Ref. Q these effects are very important as the quasi-particle thermal wave length 
is of order the average interparticle distance. 

In this paper we extend previous classical nonrelativistic simulations @ based on a color Coulomb interaction to 
the quantum regime. We develop an approach based on path integral Monte Carlo (PIMC) simulations of the strongly 
coupled QGP which self-consistently takes into account the Fermi (Bose) statistics of quarks (gluons). Following an 
idea of Kelbg quantum corrections to the pair potential are rigorously derived To extend the method of 
quantum potentials to a stronger coupling, an "improved Kelbg potential" was derived, which contains a single free 
parameter, being fitted to the exact solution of the quantum-mechanical two-body problem. Thus, the method of the 
improved Kelbg potential is able to describe thermodynamic properties up to moderate couplings [l5| . However, this 
approach may fail, if bound states of more than two particles are formed in the system. This results in a break-down 
of the pair approximation for the density matrix, as demonstrated in Ref. [l5| . A superior approach, which does not 
have this limitation, consists in use the original Kelbg potential in the PIMC simulations which effectively map the 
problem onto a high-temperature weakly coupled and weakly degenerate one. This allows one to rigorously extend 
the analysis to strong couplings and is, therefore, a relevant choice for the present purpose. 

This method has been successfully applied to strongly coupled electrodynamic plasmas [3 01- Examples are 
partially-ionized dense hydrogen plasmas, where liquid-like and crystalline behavior was observed |l8l.[l9|. Moreover, 
also partial ionization effects and pressure ionization could be studied from first principles (20l |. The same methods 
have been also applied to electron- hole plasmas in semiconductors [U, including excitonic bound states, which 
have many similarities to the QGP due to smaller mass differences as compared to electron-ion plasmas. 

The main goal of this article is to test the developed approach for ability to reproduce known lattice data d, Q and 
to predict other properties of the QGP, which are still unavailable for the lattice calculations. To this end we use a 
simple model Q of the QGP consisting of quarks, antiquarks and gluons interacting via a color Coulomb potential. 
First results of applications of the PIMC method to study of thermodynamic properties of the nonideal QGP have 
already been briefly reported in (23l [23 |. In this paper we present a comprehensive report on the thermodynamic 
properties. 

II. THERMODYNAMICS OF QGP 
A. Basics of the model 

Our model is based on a resummation technique and lattice simulations for dressed quarks, antiquarks and gluons 
interacting via the color Coulomb potential. The assumptions of the model are similar to those of Ref. Q : 

I. All color quasi-particles are heavy, i.e. their mass (m) is higher than the mean kinetic energy per particle. For 

instance, at zero net-baryon density it amounts to m > T, where T is a temperature. Therefore these particles 
move non-relativistically. This assumption is based on the analysis of lattice data [25l . [26j . 

II. Since the particles are non-relativistic, interparticle interaction is dominated by a color-electric Coulomb potential, 

see Eq. ([1]). Magnetic effects are neglected as sub- leading ones. 

III. Relying on the fact that the color representations are large, the color operators are substituted by their average 

values, i.e. by classical color vectors, the time evolution of which is described by Wong's dynamics (27| . 

The quality of these approximations and their limitations were discussed in Ref. Thus, this model requires the 
following quantities as an input: 

1. the quasi-particle mass, m, and 

2. the coupling constant g^. 

All the input quantities should be deduced from the lattice data or from an appropriate model simulating these data. 

B. Path-Integral Monte-Carlo Simulations 

Thus, we consider a three-component QGP consisting of a number of dressed quarks (Ng), antiquarks (Nq) and 
gluons (Ng) represented by quasi-particles. In thermal equilibrium the average values of these numbers can be found in 
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the grand canonical ensemble defined by the temperature-dependent Hamiltonian, which can be written as H = K+U. 
The kinetic and color Coulomb interaction energy of the quasi-particles are 



K 



l^r 1 fj l^ g^{\r^-r,lT,^i,){Q,\Q,) 



2 ^ 47r|r, 



Here the Qi denote the Wong's color variable (8-vector in the SU{i) group), T is the temperature and /ig is the quark 
chemical potential, {Qi\Qj) denote scalar product of color vectors. In fact, the quasi-particle mass and the coupling 
constant, as deduced from the lattice data, are functions of T and, in general, ^q. Moreover, is a function of distance 
r, which produces a linearly rising potential at large r (28| . 

The thermodynamic properties in the grand canonical ensemble with given temperature T, chemical potential /iq 
and fixed volume V are fully described by the grand partition function 

Z(/.„/3,F)== '^"P^^^^^ . I drdQ Pir, g, a; iV,, N,, N,; (2) 



N„,Ni,.N, 



where p{r,Q, a; Ng, Ng, Ng-, P) denotes the diagonal matrix elements of the density operator p = exp(— and 
/3 = 1/T. Here a, r and Q denote the spin, spatial and color degrees of freedom of all quarks, antiquarks and gluons 
in the ensemble, respectively. Correspondingly, the a summation and integration drdQ run over all individual degrees 
of freedom of the particles. Since the masses and the coupling constant depend on the temperature and chemical 
potential, special care should be taken to preserve thermodynamical consistency of this approach. In order to preserve 
the thermodynamical consistency, thermodynamic functions such as pressure, P, entropy, S, baryon number, Nb, and 
internal energy, E, should be calculated through respective derivatives of the logarithm of the partition function 

P = d{TlnZ)/dV, S = d{T\nZ)/dT, Ng = {l/3)d{TlnZ)/dp.g, E = -PV + TS + SfigNe- (3) 

This is a conventional way of maintaining the thermodynamical consistency in approaches of the Ginzburg-Landau 
type as they are applied in high-energy physics. 

The exact density matrix of interacting quantum systems can be constructed using a path integral approach (29l |30| 

based on the operator identity e~^^ ~ e~^^^ ■ e~^^^ . . . e~^^^ , where the r.h.s. contains n + 1 identical factors 
with A/3 = /3/(n -I- 1). which allows us to rewrite^ the integral in Eq. ^ 

J2 J dA^'^dQ^"^ p{q^°\Q^°\<j-Ng,Ng,Ng;/3)^ J dr^''^ dQ^°^ dr^^^ dQ^^K . . dA"^ dQ^'^^ p^^^ ■ p^^^ . . . 

X E E E Yi-^r^'"""'^ Sia, PgPgPgCj') PgPgPgP^-^^^ L,„+ =,,0, 
cr P, P,- Pg 

= / dQ(°)dr(0)dr(i) ...dr(")p(r(0),r(i),...r(");Q(0);7V„7V,-,7Vg;/3). (4) 



The spin gives rise to the spin part of the density matrix {S) with exchange effects accounted for by the permutation 
operators Pg, Pq and Pg acting on the quark, antiquark and gluon spatial r^""*"^) and color (5("+^) coordinates, as well as 
on the spin projections a' . The sum runs over all permutations with parity np^ and Kp-. In Eq. ^ the index 1 = 1... n-\- 

1 labels the off-diagonal density matrices = p (r^'-i) , Q('-i); r« , Q('); A/3) « (r^'^i) |e-'^'^^|r«)(5e(Q('"i) - Q(')), 
where 5,{Q^^-^^ - Q(')) is a delta-function at e — ?> 0. Accordingly each a particle is represented by a set of n -f 1 
coordinates ("beads"), i.e. by (n + 1) 3-dimensional vectors {r^\ . . . ri"^} and a 8-dimensional color vector Q*-*^^ in the 
S'C/(3) group, since all beads are characterized by the same color charge. 

The main advantage of decomposition Q is that it allows us to use a semi-classical approximation for density 
matrices which is applicable due to smallness of artificially introduced factor l/{n + 1). This parameter makes 
the thermal wavelength AAq = ■\/27rA/3/m^ of a bead of type a (a = q,q,g), smaller then a characteristic scale of 
variation of the potential energy. In the high-temperature limit p' can be approximated by a product of two-particle 
density matrices. Generalizing the electrodynamic plasma results (l7| to the case of an additional bosonic species (i.e. 



^ For the sake of notation convenience, we ascribe superscript to the original variables. 
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gluons), we write 

p(rW,rW,...r(");g(0);iV,,iV,,Ar^;/3) 
= 2^ 2N, .3N,.3N,.3N, P^^ll'/' llgluedet||0 ||.det||0 llfeHH'/'pp (5) 

where N = Ng + Ng + Ng, s and k are numbers of quarks and antiquarks, respectively, with the same 

spin projection, Aq = y^2Tr(3/ma is a thermal wavelength of an a particle, = Nal/[s\{Na — s)!], 

the antisymmetrization and symmetrization are taken into account by the symbols "det" and "per" denoting 



the determinant and permanent, respectively. Functions <j>pp = exp 



and matrix elements 0"^*^ 



2 



_ ^,{0)-^ y^"^ /AA^ ) S^{Qt — Qo), where t and o are particle's indexes, are expressed in terms of 



exp 

distances . . . , and dimensionless distances {£,i^\ ■ ■ ■ , ) between neighboring beads of an a particle, 

defined as Tq'' = ri"' + ya\ (/ > 0), and ya"* = AAq ^[_]^ ^i*'"'. Notice that the indices s and k in det||0"'°||s 
and det ||0"'°||k denoted that matrices ||0"'°|| have two nonzero blocks related to quark and antiquark quasi-particles 
with the same spin projections. The density matrix ([S]) has been transformed to a form which does not contain an 
explicit sum over permutations and thus no sum of terms with alternating sign (in the case of quarks and antiquarks) . 
Let us stress that the determinants depend also on the color variables. 
In Eq. ^ the total color interaction energy 

-, n + l -, n+1 -| 

u{r, Q,p) = -L-^Y. u''' = E ^ E - -t"^ I ' l'-^'^ - ^t^lQp^ Q^) (6) 

1=1 1=1 p^t 

is defined in terms of off-diagonal two-particle effective quantum potential which is obtained by expanding the 
two-particle density matrix ppt up to the first order in small parameter l/{n + 1): 

Jo J 47r|r"|AA2j^r(f -t) 

"""P (~ AAg,(f -''') ) """P (' '^AA^.r' ) P°*exp[-A/3<i>^'*(r,r',g„Q0]- (7) 



where r = rp — rt, i^' = r'p ^ r'^, AApt — ^27rA/3/TOp4, uipt = mpnit / (mp + rrit) is the reduced mass of the (j)i)-pair of 
particles, and p^j is the two-particle density matrix of the ideal gas. The result for the diagonal color Kelbg potential 
can be written as 

r, Qp, Q,) « 9\T^,\{Qp\Q,) r _ ^ _ )j\ 

ATT/\\ptXpt L J 



where Xpt = \rp — rt|/AApt. Here the function g (T,fiq) = g^(r", T./i^), resulting from averaging of the initial 
g'^{r" , T, jig) over relevant distances of order AApt, plays the role of an effective coupling constant. Note that the color 
Kelbg potential approaches the color Coulomb potential at distances larger than A Apt. What is of prime importance, 
the color Kelbg potential is finite at zero distance, thus removing in a natural way the classical divergences and making 
any artificial cut-offs obsolete. This potential is a straightforward generalization of the corresponding potential of 
electrodynamic plasmas [T5|. The off-diagonal elements of the effective interaction are approximated by the diagonal 
one by means of $P*(r, r'; , Qp, Qt) « [^P\r,r,Qp,Qt) + <^P\r' y ,Qp,Qt)]/2. 

The described path-integral representation of the density matrix is exact in the limit n — > oo. For any finite number 
n, the error of the above approximations for the whole product on the r.h.s. of Eq. @ is of the order f /(?! + f ) whereas 
the error of each p' is of the order f /(n -t- 1)^, as it was shown in Ref. [T7| . 

The main contribution to the partition function comes from configurations in which the 'size' of the cloud of beads 
of quasi-particles is of the order of their thermal wavelength Aq whereas characteristic distances between beads of 
each quasi-particle are of the order of AAq. 



III. SIMULATIONS OF QGP 

To test the developed approach we consider the QGP only at zero baryon density and further simplify the model 
by additional approximations, similarly to Ref. Q: 
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IV We replace the grand canonical ensemble by a canonical one. The thermodynamic properties in the canonical 

ensemble with given temperature T and fixed volume V are fully described by the density operator p = e~^^ 
with the partition function defined as follows 

Z{N,,N,,Ng,V;(3) = ^ , / drdQ p{r,Q,a; f3), (9) 

9 « 9 o- ^ 

with Nq = Nq and hence Nb = 0. In order to preserve the thermodynamical consistency of this formulation, 
thermodynamic quantities should be calculated through respective derivatives of the logarithm of the partition 
function similarly to that in Eq. ^ with the exception that now Na are independent variables. 

V Since the masses of quarks of different flavors extracted from lattice data are very similar, we do not distinguish 

between quark flavors. Moreover, we take the quark and gluon quasi-particle masses being equal because their 
values deduced from the lattice data (25l . [26| are very close. 

VI Because of the equality of masses and approximate equality of number of degrees of freedom of quarks, antiquarks 

and gluons, we assume that these species are equally represented in the system: Nq = Nq ~ Ng. 

VII For the sake of technical simplicity, the SU{3) color group is replaced by SU{2). 
Thus, this simplified model requires an additional quantity as an input: 

3. the density of quasi-particles (Nq + Ng + Ng)/V = n{T) as a function of the temperature. 

Although this density is unknown from the QCD lattice calculations and we use it as a fit parameter, it is very 
important to partially overcome constrains of the above simplifications. First, it concerns the use of the SU{2) color 
group, which first of all reduces the degeneracy factors of the quark and gluon states, as compared to the SU{3) case, 
and thereby reduces pressure and all other thermodynamic quantities. A proper fit of the density allows us to remedy 
this deficiency of the normalization. Second, in fact we consider the system of single quark flavor, i.e. all quarks are 
identical, which also reduces the normalization of all thermodynamic quantities. The density fit cures the deficiency 
of this normalization, though the excessive anticorrelation of quarks remains. 

Ideally the parameters of the model should be deduced from the QCD lattice data. However, presently this task is 
still quite ambiguous. Thereforejin the present simulations we take a possible (maybe, not the most reliable) set of 
parameters. Following Refs. (ol. [26|. the parametrization of the quasi-particle mass is taken in the form 

m(T) /Tc = 0.9 /{T/Tc - 1) + 3.45 + OAT/Tc (10) 

where Tc = 175 MeV is the critical temperature. This parametrization fits the quark mass at two values of temperature 
obtained in the lattice calculations [1^. According to [2^ the masses are quite large: niq/T ~ 4 and rag/T ~ 3.5. These 
are essentially larger than masses required for quasi-particle fits (3ll |32| of the lattice thermodynamic properties of 
the QGP: niq/T ~ 1-^2 and irig/T ~ 1.5-;- 3. Moreover, the pole quark mass rUq/T ~ 0.8 was reported in recent work 
(33j . as deduced from lattice calculations. Nevertheless, in spite of the fact that it obviously produces too high masses, 
we use parametrization (jTU]) in order to be compatible with the input of classical MD of Ref. Q. The T-dependence 
of this mass is illustrated in Fig. [T] (left panel). 

The coupling constant, i.e. as = /{Att), used in the simulations is displayed in the left panel of Fig. [T] as well. As 
seen, cts well complies with phenomenologic QCD estimations [34l] of its values. Notice that in previous publications 
(23l. the factor {N'^ — 1), where Nc is the number of colors in the SU{Nc) color group, was accidentally included 
in g^, when displaying it corresponding figures. In fact, this factor is a part of Casimirs defining the normalization of 
the color vectors in the color group. 

The density of quasi-particles, which is additionally required within the canonical-ensemble approach, was chosen 
on the condition of the best agreement of the calculated pressure with the corresponding lattice result, see Fig. [T] 
(right panel). It was taken to be n{T) = 0.24r'^. From the first glance, it is a very low density. For example, in the 
classical simulations of Ref. Q it was taken as n{T)/T^ = 6.3, which corresponds to the density of an ideal gas of 
massless quarks, antiquarks and gluons. Since the quasi-particles are very heavy in the present model (as well as in 
that of Ref. [^), the latter density looks unrealistically high. Even in quasi-particle models [Sllll^l, where the masses 
are lower, the density turns out to be n{T)/T^ « 1.4. Since Eq. ([TU)) gives even larger masses than those in Refs. 
(31I . [3^ and in view of the adopted large coupling, the chosen value of n{T) does not look too unrealistic. 

Thus, although the chosen set of parameters is still debatable, it is somehow self-consistent. In the future we are 
going to get rid of the n{T) parameter, by applying the grand-canonical approach, and by using more moderate (and 
maybe realistic) sets of parameters. 
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Calculation of the equation of state (right panel of Fig. [T]) was used to optimize the parameters of the model in order 
to proceed to predictions of other properties concerning the internal structure and in the future also non-equilibrium 
dynamics of the QGP. The plasma coupling parameter is defined as 

r = ^ (11) 

where Ts is the the Wigner-Seitz radius, defined such that 47rr^/3 — n, and 02 the quadratic Casimir value averaged 
over quarks, antiquarks and gluons. 02 = — 1 is a good estimate for this quantity. The plasma parameter is a 
measure of ratio of the average potential to the average kinetic energy. It is also presented in Fig. [1] (left panel). It 
turns out to be of the order of unity which indicates that the QGP is a strongly coupled Coulomb liquid rather than a 
gas. In the studied temperature range, 1 < T/Tc < 3, the QGP is, in fact, quantum degenerate, since the degeneracy 
parameter Xa = "-aAj^ (where the thermal wave length. A, is defined in the previous sect.) varies from 0.1 to 1.7, see 
Fig. □ (left panel). 




Phc. 1: Left panel: Temperature dependence of the model input quantities, coupling constant and mass-to-temperature ratio 
(scaled by 1/2), the plasma coupling parameter F [see Eq. (|lip ] and the degeneracy parameter x- The x parameters for different 
species are equal, since their masses and densities are assumed to be equal. Right panel: Equation of state (pressure versus 
temperature) of the QGP from PIMC simulations (open squares) compared to lattice data of Refs. 0, The solid line is a 
smooth interpolation between the PIMC points. Results of the HotQCD Collaboration Q are presented by filled circles, while 
results of the Budapest group open circles. Different kinds of circles (filled and open) correspond to different discretization 
schemes of the QCD action (p4 and asqtad, see [3| for details). 




Phc. 2: Entropy (left panel) and trace anomaly (right panel) of the QGP from PIMC simulations (solid line) compared to 
lattice data of Refs. Q, Q- Notation is the same as in the right panel of Fig. [T] 

Figure [2] additionally presents the entropy (S) and trace anomaly (e — 3P) of the QGP computed in the PIMC 
method. These quantities are calculated accordingly to Eqs. In order to avoid the numeric noise, the derivative of 
a smooth interpolation between the PIMC points (see the right panel of Fig.[T]) was taken. These results are compared 
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to lattice data of Refs. 0, It is not surprising that agreement with the lattice data is also good, since it is a direct 
consequence of the good reproduction of the pressure. 

Details of our pat h integral Monte-Carlo simulations have been discussed elsewhere in a variety of papers and review 
articles, see, e.g. [33 and references therein. The main idea of the simulations consists in constructing a Markov process 
of configurations which differ by the particle coordinates. Additionally to the case of electrodynamic plasmas, here 
we randomly sample, according to the group measure, the color variables Q of all particles until convergence is 
achieved. We use a cubic simulation box with periodic boundary conditions. The number of particles was taken as 
N = Nq + Nq + Ng = + + AO ^ 120, and the number of beads, n = 20. Calculation of the equation of state 
(right panel of Fig. [T]) was used to optimize the parameters of the model in order to proceed to predictions of other 
properties concerning the internal structure and in the future also non-equilibrium dynamics of the QGP. 





Phc. 3: Pair distribution functions (upper panels) and color pair distribution functions (lower panels) of identical (left panels) 
and different (right panels) quasi-particles at temperature T /Tc — 3. The distance is measured in units of ct = 1/Tc = 1.1 fm. 

Let us now consider the spatial arrangement of the quasi-particles in the QGP by studying the pair distribution 
functions (PDF's) gabif)- They give the probability density to find a pair of particles of types a and 6 at a certain 
distance r from each other and are defined as 

9abiRi - R2) = ,^ ,^ , I ^'^^Q ^(^1 - - r^)p(r, g, a; (3). (12) 

Q Q 3 ^ 

The PDF's depend only on the difference of coordinates because of the translational invariance of the system. In a 
non-interacting classical system, gab{r) = 1, whereas interactions and quantum statistics result in a re-distribution of 
the particles. Results for the PDF's at temperature T/Tc ~ 3 are shown in the top panels of Fig. [3l 

At large distances, r/a > 0.5, all PDF's of identical particles (top left panel of Fig.[3|) coincide. A drastic difference 
in the behavior of the PDF's of quarks and gluons (the anti-quark PDF is identical to the quark PDF) occurs at 
small distances. Here the gluon PDF increases monotonically when the distance goes to zero, while the PDF of quarks 
(and antiquarks) exhibits a broad maximum. This difference is the effect of the quantum statistics. In the present 
conditions, the thermal wavelength A approximately equals 0.37cr, i.e. the difference starts to appear at distances 
of the order of A. The enhanced population of low distance states of gluons is due to Bose statistics and the color- 
Coulomb attraction. In contrast, the depletion of the small distance range for quarks is a consequence of the Pauli 
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principle. In an ideal Fermi gas g{r) equals zero for particles with the same spin projections and colors, while for 
particles with different colors and/or opposite spins the PDF equals unity in the limit r — >■ 0. As a consequence, 
the spin and color averaged PDF approaches 0.5. Such a low-distance behavior is also observed in a nonideal dense 
astrophysical electron-ion plasma and in a nonideal electron- hole plasmas in semiconductors (2ll . [2^ . The depletion 
of the probability of quasi-particles at small distances results in its enhancement at intermediate distances. This is 
the reason for the corresponding PDF maxima. Oscillations of the PDF at very small distances of order r < 0.02(t are 
related to Monte Carlo statistical error, as probability of quasi-particles being at short distances quickly decreases. 

At small distances, r < 0.3cr, a strong increase is observed in all PDF's of particles of different type (top right panel 
of Fig. [3]), which resembles the behavior of the gluon-gluon PDF. This increase is a clear manifestation of an effective 
pair attraction of quarks and antiquarks as well as quarks (antiquarks) and gluons. This attraction suggests that the 
color vectors of nearest neighbors of any type are anti-parallel. If this explanation is correct can be verified by means 
of color pair distribution functions (CPDF) defined as 

Cab{Ri - i?2) = ,^ ,^ , E / ^''^Q mQ\)mi - r1)SiR2 - rl)p{r, Q, a; /3), (13) 

which are shown in the lower panels of Fig. [31 All CPDF's turn out to be negative at small distances, indicating 
anti-parallel orientation of the color vectors of neighboring quasi-particles, which complies with attraction seen in the 
functions gab for a ^ b. The minimum of Cqq close to r = 0.2a corresponds to the maximum observed in gqq. The deep 
minimum in the gluon CPDF at small distances results from the Bose statistics and complies with the high maximum 
of the gluon PDF ggg. 

Thus, at T/Tc = 3 we observe signs of a spatial ordering, cf. the peak of the quark PDF around r/cr = 0.1 — 0.2, 
which may be interpreted as emergence of liquid-like behavior of the QGP. Much more pronounced is the short-range 
structure of nearest neighbors. The QGP lowers its total energy by minimizing the color Coulomb interaction energy 
via a spontaneous "anti-ferromagnetic" ordering of color vectors. This gives rise to a clustering of quarks, antiquarks 
and gluons. To verify the relevance of these trends, a more refined spin-resolved analysis of the PDF's and CPDF's is 
necessary, together with simulations in a broader range of temperatures which are presently in progress. 

PDF's allow us also to do a lot of physical conclusions rather than only to analyze the spatial structure of the 
QGP. Fig. |4] presents PDF's of the identical particles for two temperatures T = l.lTc and T = 2Tc (left upper and 
middle panels). These PDF's can be formed either by correlated scattering states or by bound states of quasi-particles, 
depending on the relative fractions of these states. In a strict sense, however, there is no clear subdivision into bound 
and free "components" due to the mutual overlap of the quasi-particle clouds. In addition, there exists no rigorous 
criterion for a bound state at high densities due to the strong effect of the surrounding plasma. Nevertheless, a rough 
estimate of the fraction of quasi-particle bound states can be obtained by the following reasonings. The product 
r^gabir) has the meaning of a probability to find a pair of quasiparticles at a distance r from each other. On the other 
hand, the corresponding quantum mechanical probability is the product of and the two-particle Slater sum 

l^ab = Sir'^'-Xlb E I*" Wl' ^M-PEa.) - ^ib + Kb, (14) 

a 

where Ea and 5'a(r) are the energy (without center of mass energy) and the wave function of a quasi-particle pair. , 
respectively, and Xab = \/27r/3(ma + rub) / {mamb). Sab is, in essence, the diagonal part of the corresponding density 
matrix. The summation runs over all states a of the discrete (SJ^b) ^^^"^ continuous (Kb) spectrum. The discrete- 
spectrum contribution 

E' 

= 8^3/2 A^, ^ |vI/„(r)pexp(-/?£;„) (15) 

E^=Eo 

originates from the populated states between the ground one with the energy Eq and and some state with the energy 
E' . For low densities it is reasonable to choose E' sa —T (30j . 

At temperatures smaller than the binding energy and distances smaller than or of the order of several bound state 
radii, the main contribution to the Slater sum comes from bound states, while at larger distances free (scattering) 
states give a dominant contribution. In the electromagnetic plasma it was found that the product r'^'Sff^ is sharply 
peaked at distances around the Bohr radius in this case. Similarly, at low temperature, r^gqq(r) forms a pronounced 
maximum near r = 0.2 fm which can be interpreted as the radius of a bound qq pair (see right upper panel of Fig. 
HI). Thus, our calculations support the existence of bound states of medium- modified (massive) quarks and gluons at 
moderate temperatures, i.e. just above T^, proposed in Ref. [sBl and later in Refs. [s^jllll based on results from lattice 
QCD calculations of spectral functions |39l |40|. Integration of Cqq from zero up to 0.3 fm shows that modulus scalar 
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Phc. 4: Top and middle panels: Pair distribution functions at two different temperatures T (left panels) and quark-antiquark 
PDF multiplied by distance squared (right panels). The distance is measured in units of cr = 1/Ta = 1.1 fm. Bottom panels: 
potential of average force, defined according to Eq. (|16|) . at small distances for T — l.lTc (right panel) and logarithm of the 
product r and averaged potential at large distances (left panel). Dashed lines are linear fits to the calculated (solid) lines. 

product \{Qq + QqlQq + Qq)\ IS much smaller that \{Qq\Qq)\ and that means that these bound states are 

colorless. With the temperature rise these bound states dissolve much faster than it was assumed in fs?!. Issj. which 
complies with the analysis of Ref. j4l| . Indeed, at the temperature of T = 2Tc the bound states completely disappear 
(see right middle panel of Fig. |4]) and r^gatif) approaches the behaviour which corresponds to the plasma without 
bound states. 

Interesting observations can be done from the analysis of the potential of average force (PAF) defined as the 
logarithm of the related PDF, 



Uab{r,T) = ~T\ngab{r,T) 



(16) 
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(lower panels of Fig. H]). This definition is motivated by the PDF virial expansion in terms of bare potential (like color 
Kelbg potential) . Near the QGP phase transition the PAF (right bottom panel) is a linear function at distances smaller 
than the bound state radius. This suggests that the bound states are bound by a string-like forces. At larger distances 
the PAF (left bottom panel) can be very well approximated by an exponentially screened Coulomb potential (Yukawa- 
type potential) like that in the electromagnetic plasma. Potential well of the linear part of average force potential is 
of order 1.5 GeV at T = l.lTc, while for T = 2Tc it is only of order 0.2 Gev. Straight dashed lines in the lower panels 
of Fig. m are the least-square liner (left panel) and Yukawa-type (right panel) approximations to the the effective 
potential. 



IV. CONCLUSION 



Quantum Monte Carlo simulations based on the quasi-particle picture of the QGP are able to reproduce the lattice 
equation of state (even near the critical temperature) and also yield valuable insight into the internal structure of 
the QGP. Our results indicate that the QGP reveals liquid-like (rather than gas-like) properties even at the highest 
considered temperature of STc. At temperatures just above Tc we have found that bound quark- ant iquark states still 
survive. These states are bound by effective string-like forces. Quantum effects turned out to be of prime importance 
in these simulations. 

Our analysis is still too simplified and incomplete. It is still confined only to the case of zero baryon chemical 
potential. The input of the model also requires refinement. Work on these problems is in progress. 

However, the PIMC method is not able to yield dynamical and transport properties of the QGP. One way to 
achieve this is to develop semiclassical color molecular dynamics simulations. In contrast to previous MD simulations 
0, where quantum effects were included phenomenologically via a short range potential, we have developed a more 
rigorous approach to study the dynamical and transport properties of strongly coupled quark-gluon systems which 
is based on the combination of the Feynman and Wigner formulation of quantum dynamics. The basic ideas of this 
approach for the electron- ion plasmas have been published in [T^ . [Tsl . [35| . In particular, this approach allows us to 
deduce the viscosity of the QGP. Work on these problems is also in progress. 

We acknowledge stimulating discussions with D. Blaschke, M. Bleicher, R. Bock, B. Friman, C. Ewerz, D. Rischke, 
and H. Stoecker. Y.I. was partially supported by the Deutsche Forschungsgemeinschaft (DFG projects 436 RUS 
113/558/0-3 and WA 431/8-1), the RFBR grant 09-02-91331 NNIO_a, and grant NS-7235.2010.2. 
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